%     VIRTUAL FABRIC - программа для получения эффективных свойств текстильных композитов
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

clear all; %#ok<CLALL>
clc;
more off;

% #######################################################################################################
% ##                                                                                                   ##
% ##                                        НАСТРОЙКИ ПРОГРАММЫ                                        ##
% ##                                                                                                   ##
% #######################################################################################################

% Для включения возможности вывода информации на русском языке при выполнении программы в среде GNU Octave
% необходимо добавить в файл C:\Octave\Octave-4.4.0\share\octave\site\m\startup\octaverc  следующую строку:
% system("chcp 1251",1);

% --------------------------------- Параметры воксельной модели TexGen ----------------------------------

Voxel_Model_Directory      = 'M:\Abaqus_Work';          % Папка с воксельной моделью, созданной в TexGen
Voxel_Model_Name           = 'fabric_voxel_example';    % Название воксельной модели (имя .inp файла)

% --------------------------------- Параметры объемной доли связующего ----------------------------------

V_resin_Type               = 'LOCAL'; % LOCAL | GLOBAL    % Тип объемной доли связующего: LOCAL - локальная объемная доля в жгутах, GLOBAL - суммарная объемная доля в ткани.
V_resin                    = 0.25;    % Объемная доля связующего

% ----------------------------------------- Параметры расчета -------------------------------------------

Analysis_Type              = 'MECHANICAL'; % MECHANICAL | HEAT_TRANSFER | MESH_GENERATOR

                           % MECHANICAL - получение эффективных механических свойств с помощью решателя Abaqus
                           % HEAT_TRANSFER - получение эффективных свойств теплопереноса с помощью решателя Abaqus
                           % MESH_GENERATOR - только генерация сетки (без расчета)
                           
Mechanical_Element_Type    = 'C3D8I'; % C3D8R | C3D8RH | C3D8 | C3D8H | C3D8I | C3D8IH    % Тип конечного элемента в воксельной модели для расчета механических свойств
                           
Mechanical_BC_TYPE         = 'SOFT'; % SOFT | HARD    % Тип граничных условий в механическом расчете

                           % SOFT - периодические граничные условия со свободной границей ячейки
                           % HARD - периодические граничные условия с фиксированной границей ячейки (использовать с большой осторожностью!)
                           
Constraint_Directions      = 'XYZ';    % Направления, задающие грани с фиксированной границей в случае Mechanical_BC_TYPE == HARD. Строка, состоящая из символов X, Y, Z.

Chemical_Shrinkage         = 'FALSE'; % TRUE | FALSE    % Расчет эффективных коэффициентов химической усадки помимо основных механических свойств

% ----------------------------------------- Свойства материала ------------------------------------------

% ДЛЯ МЕХАНИЧЕСКОГО И ТЕПЛОВОГО РАСЧЕТА

Fiber_rho                  = 1.78;    % Плотность волокна
Matrix_rho                 = 1.2;    % Плотность матрицы

% ДЛЯ МЕХАНИЧЕСКОГО РАСЧЕТА

% Свойства волокна
Fiber_E_L                  = 240;    % Модуль упругости в продольном направлении (L)
Fiber_E_T                  = 15;    % Модуль упругости в поперечном направлении (T)
Fiber_G_LT                 = 26;    % Модуль сдвига в плоскости LT
Fiber_nu_LT                = 0.267;    % Коэффициент Пуассона LT (nu_LT*E_T = nu_TL*E_L)
Fiber_nu_TT                = 0.485;    % Коэффициент Пуассона TT
Fiber_alpha_L              = -0.568;    % Коэффициент линейного теплового расширения в направлении L
Fiber_alpha_T              = 9.372;    % Коэффициент линейного теплового расширения в направлении T
Fiber_beta_L               = 0;    % Коэффициент химической усадки в направлении 1
Fiber_beta_T               = 0;    % Коэффициент химической усадки в направлении 2

% Свойства матрицы
Matrix_E                   = 3.5;    % Модуль упругости
Matrix_nu                  = 0.38;    % Коэффициент Пуассона
Matrix_alpha               = 70;    % Коэффициент линейного теплового расширения
Matrix_beta                = -0.01;    % Коэффициент химической усадки

% ДЛЯ ТЕПЛОВОГО РАСЧЕТА

% Свойства волокна
Fiber_Cp                   = 710;    % Удельная теплоемкость
Fiber_lambda_L             = 10;    % Коэффициент теплопроводности в продольном направлении (L)
Fiber_lambda_T             = 1;    % Коэффициент теплопроводности в поперечном направлении (T)

% Свойства матрицы
Matrix_Cp                  = 1200;    % Удельная теплоемкость
Matrix_lambda              = 0.2;    % Коэффициент теплопроводности

% ----------------------------------------- Параметры запуска -------------------------------------------

Work_Directory             = 'M:\Abaqus_Work';    % Рабочая папка для проведения расчета в Abaqus
Solver_Command             = 'abaqus';            % Команда вызова решателя. Примеры: abaqus, abq2016hf7, abq2019hf9
CPU                        = 4;                   % Суммарное число используемых процессорных ядер для расчета в Abaqus
DMP                        = 1;                   % Число локальных DMP-процессов (по CPU/DMP потоков в каждом)
Memory                     = [];                  % Максимальный объем оперативной памяти, используемой в расчете воксельной модели. Пример: 'memory="220 Gb"'

Mech_Solver                = 'ITERATIVE';         % DIRECT | ITERATIVE    % Тип решателя в расчете механических свойств на воксельной модели
Mech_Max_Iter              = 500;                 % Максимальное число итераций (только для итерационного решателя)
Mech_Rel_Tol               = 1e-6;                % Величина относительной невязки (только для итерационного решателя)

Heat_Transfer_Solver       = 'DIRECT';            % DIRECT | ITERATIVE    % Тип решателя в расчете свойств теплопереноса на воксельной модели
Heat_Transfer_Max_Iter     = 500;                 % Максимальное число итераций (только для итерационного решателя)
Heat_Transfer_Rel_Tol      = 1e-6;                % Величина относительной невязки (только для итерационного решателя)

% -------------------------------------------------------------------------------------------------------

% #######################################################################################################
% ##                                                                                                   ##
% ##                                     КОНЕЦ НАСТРОЕК ПРОГРАММЫ                                      ##
% ##                                                                                                   ##
% #######################################################################################################

global ErrCode;
global Merge_Tolerance_Multiplier;
global software_version;
global software_date;
global software_year;

global Auto_Thickness_Multiplier;
global Mesh_Seed_M_Multiplier;

ErrCode = 0;
Merge_Tolerance_Multiplier = [];
software_version = [];
software_date = [];
software_year = [];

Auto_Thickness_Multiplier = [];
Mesh_Seed_M_Multiplier = [];

% Отсчет времени выполнения программы
tic;

% Инициализация базовых переменных
Init_Variables();

% Печать приветствия
Print_Greeting();

% Проверка параметров
switch V_resin_Type
    case {'LOCAL', 'GLOBAL'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра V_resin_Type !\n');
        fprintf (1, 'Разрешенные значения: LOCAL, GLOBAL\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if ((V_resin <= 1 - pi()/(2*sqrt(3))) || (V_resin >= 1))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра V_resin !\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Analysis_Type
    case {'MECHANICAL', 'HEAT_TRANSFER', 'MESH_GENERATOR'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Analysis_Type !\n');
        fprintf (1, 'Разрешенные значения: MECHANICAL, HEAT_TRANSFER, MESH_GENERATOR\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

switch Mechanical_Element_Type
    case {'C3D8R', 'C3D8RH', 'C3D8', 'C3D8H', 'C3D8I', 'C3D8IH'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Mechanical_Element_Type !\n');
        fprintf (1, 'Разрешенные значения: C3D8R, C3D8RH, C3D8, C3D8H, C3D8I, C3D8IH\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

switch Mechanical_BC_TYPE
    case {'SOFT', 'HARD'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Mechanical_BC_TYPE !\n');
        fprintf (1, 'Разрешенные значения: SOFT, HARD\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if ((numel(Constraint_Directions) < 1) || (numel(Constraint_Directions) > 3))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Constraint_Directions !\n');
    fprintf (1, 'Разрешенные значения: строка, состоящая из символов X, Y, Z.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

for i=1:numel(Constraint_Directions)
    switch Constraint_Directions(i)
        case {'X', 'Y', 'Z'}
            % Допустимые значения
        otherwise
            ErrCode = -1;
            fprintf (1, 'Недопустимое значение параметра Constraint_Directions !\n');
            fprintf (1, 'Разрешенные значения: строка, состоящая из символов X, Y, Z.\n');
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
    end
end

switch Chemical_Shrinkage
    case {'TRUE', 'FALSE'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Chemical_Shrinkage !\n');
        fprintf (1, 'Разрешенные значения: TRUE, FALSE\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (CPU ~= floor(CPU))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра CPU !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (CPU < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра CPU !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (DMP ~= floor(DMP))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра DMP !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (DMP < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра DMP !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Mech_Solver
    case {'DIRECT', 'ITERATIVE'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Mech_Solver !\n');
        fprintf (1, 'Разрешенные значения: DIRECT, ITERATIVE\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (Mech_Max_Iter ~= floor(Mech_Max_Iter))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Mech_Max_Iter !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (Mech_Max_Iter < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Mech_Max_Iter !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (Mech_Rel_Tol <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Mech_Rel_Tol !\n');
    fprintf (1, 'Параметр должен быть больше нуля.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

switch Heat_Transfer_Solver
    case {'DIRECT', 'ITERATIVE'}
        % Допустимые значения
    otherwise
        ErrCode = -1;
        fprintf (1, 'Недопустимое значение параметра Heat_Transfer_Solver !\n');
        fprintf (1, 'Разрешенные значения: DIRECT, ITERATIVE\n');
        fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
        toc;
        return;
end

if (Heat_Transfer_Max_Iter ~= floor(Heat_Transfer_Max_Iter))
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Heat_Transfer_Max_Iter !\n');
    fprintf (1, 'Параметр должен быть целым числом.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (Heat_Transfer_Max_Iter < 1)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Heat_Transfer_Max_Iter !\n');
    fprintf (1, 'Параметр должен быть больше или равен единице.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

if (Heat_Transfer_Rel_Tol <= 0)
    ErrCode = -1;
    fprintf (1, 'Недопустимое значение параметра Heat_Transfer_Rel_Tol !\n');
    fprintf (1, 'Параметр должен быть больше нуля.\n');
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Счетчики числа сообщений
N_Info = 0;
N_Warnings = 0;

% Анализ настроек программы
if (strcmp(Analysis_Type, 'MECHANICAL'))
    
    if (strcmp(Mechanical_BC_TYPE, 'HARD'))
        fprintf (1, 'ПРЕДУПРЕЖДЕНИЕ! В расчете используются периодические граничные условия с фиксированной границей\n');
        fprintf (1, '                ячейки (Mechanical_BC_TYPE = HARD). Требуется особое внимание при анализе\n');
        fprintf (1, '                полученных эффективных свойств.\n\n');
        N_Warnings = N_Warnings + 1;
    end
    
end

% Корректировка Voxel_Model_Directory
if (Voxel_Model_Directory(length(Voxel_Model_Directory)) == '\')
    Voxel_Model_Directory = Voxel_Model_Directory(1:length(Voxel_Model_Directory)-1);
end

% Корректировка Work_Directory
if (Work_Directory(length(Work_Directory)) == '\')
    Work_Directory = Work_Directory(1:length(Work_Directory)-1);
end

% Корректировка имени inp файла
if ((length(Voxel_Model_Name) > 4) && strcmpi(Voxel_Model_Name(length(Voxel_Model_Name)-3:length(Voxel_Model_Name)), '.inp'))
    Voxel_Model_Name = Voxel_Model_Name(1:length(Voxel_Model_Name)-4);
end

% Заполнение структуры со свойствами материала
Material_Properties.Fiber_E_1 = Fiber_E_L;
Material_Properties.Fiber_E_2 = Fiber_E_T;
Material_Properties.Fiber_E_3 = Fiber_E_T;
Material_Properties.Fiber_G_12 = Fiber_G_LT;
Material_Properties.Fiber_G_13 = Fiber_G_LT;
Material_Properties.Fiber_G_23 = Fiber_E_T/2/(1+Fiber_nu_TT);
Material_Properties.Fiber_Nu_12 = Fiber_nu_LT;
Material_Properties.Fiber_Nu_13 = Fiber_nu_LT;
Material_Properties.Fiber_Nu_23 = Fiber_nu_TT;
Material_Properties.Fiber_alpha_11 = Fiber_alpha_L;
Material_Properties.Fiber_alpha_22 = Fiber_alpha_T;
Material_Properties.Fiber_alpha_33 = Fiber_alpha_T;
Material_Properties.Fiber_beta_11 = Fiber_beta_L;
Material_Properties.Fiber_beta_22 = Fiber_beta_T;
Material_Properties.Fiber_beta_33 = Fiber_beta_T;
Material_Properties.Fiber_rho = Fiber_rho;
Material_Properties.Fiber_Cp = Fiber_Cp;
Material_Properties.Fiber_lambda_11 = Fiber_lambda_L;
Material_Properties.Fiber_lambda_22 = Fiber_lambda_T;
Material_Properties.Fiber_lambda_33 = Fiber_lambda_T;

Material_Properties.Matrix_E = Matrix_E;
Material_Properties.Matrix_nu = Matrix_nu;
Material_Properties.Matrix_alpha = Matrix_alpha;
Material_Properties.Matrix_beta = Matrix_beta;
Material_Properties.Matrix_rho = Matrix_rho;
Material_Properties.Matrix_Cp = Matrix_Cp;
Material_Properties.Matrix_lambda = Matrix_lambda;

% Проверка свойств материала
Check_Material_Properties(Material_Properties);

if (ErrCode < 0)
    fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
    toc;
    return;
end

% Заполнение структуры с параметрами решателя
Solver_Parameters.Mech_Solver = Mech_Solver;
Solver_Parameters.Mech_Max_Iter = Mech_Max_Iter;
Solver_Parameters.Mech_Rel_Tol = Mech_Rel_Tol;

Solver_Parameters.Heat_Transfer_Solver = Heat_Transfer_Solver;
Solver_Parameters.Heat_Transfer_Max_Iter = Heat_Transfer_Max_Iter;
Solver_Parameters.Heat_Transfer_Rel_Tol = Heat_Transfer_Rel_Tol;

% Ветвление по типу анализа
switch Analysis_Type
    
    case 'MECHANICAL'
        
        Mechanical_Procedure(Voxel_Model_Directory, Voxel_Model_Name, Analysis_Type, Mechanical_BC_TYPE, Constraint_Directions, Chemical_Shrinkage, ...
                             Work_Directory, CPU, DMP, Material_Properties, V_resin_Type, V_resin, Memory, Solver_Command, Solver_Parameters, Mechanical_Element_Type);
        
        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
        
    case 'HEAT_TRANSFER'
        
        Heat_Transfer_Procedure(Voxel_Model_Directory, Voxel_Model_Name, Analysis_Type, Work_Directory, CPU, DMP, Material_Properties, V_resin_Type, ...
                                V_resin, Memory, Solver_Command, Solver_Parameters, Mechanical_Element_Type);
        
        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
        
    case 'MESH_GENERATOR'
        
        Mesh_Generator_Procedure(Voxel_Model_Directory, Voxel_Model_Name, Analysis_Type, Work_Directory, Material_Properties, ...
                                 V_resin_Type, V_resin, Mechanical_Element_Type);

        if (ErrCode < 0)
            fprintf (1, 'Работа программы завершена с ошибкой!\n\n');
            toc;
            return;
        end
    
end

% Печать разделителя
Print_Delimiter();

% Сообщение об успешном завершении работы программы
fprintf (1, 'Работа программы завершена успешно!\n\n');

% Вывод числа предупреждений
if (N_Warnings > 0)
    fprintf (1, '\tПолучено предупреждений: %i\n', N_Warnings);
end

% Вывод числа информационных сообщений
if (N_Info > 0)
    fprintf (1, '\tПолучено информационных сообщений: %i\n', N_Info);
end

% Печать новой строки
if ((N_Warnings > 0) || (N_Info > 0))
    fprintf (1, '\n');
end

% Вывод времени выполнения программы
toc;
